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Ultrasound is an important imaging modality for biological systems. High-frequency ultrasound 
can also (e.g., via acoustical nonlinearities) be used to provide deeply penetrating and high-resolution 
imaging of vascular structure via catheterisation. The latter is an important diagnostic in vascular 
health. Typically, ultrasound requires sources and transducers that are greater than, or of order 
the same size as the wavelength of the acoustic wave. Here we design and theoretically demonstrate 
that single silver nanorods, acting as optical nonlinear dipole antennae, can be used to detect 
ultrasound via Brillouin light scattering from linear and nonlinear acoustic waves propagating in 
bulk water. The nanorods are tuned to operate on high-order plasmon modes in contrast to the 
usual approach of using fundamental plasmon resonances. The high-order operation also gives rise 
to enhanced optical third-harmonic generation, which provides an important method for exciting 
the higher-order Fabry-Perot modes of the dipole antenna. 


I. INTRODUCTION 

Control of light at the nanoscale is a burgeoning held, 
and optical antennae are fast emerging as one of the most 
important tools for collecting, emitting and more gen¬ 
erally controlling light with nanoscale (sub wavelength) 
elements Qi- This functionality derives from strong 
coupling of optical helds to plasmon modes in the anten¬ 
nae, where the resonance of the plasmon modes can be 
controlled by the geometry of the nanoscale elements. 

We are interested in the detection of high-frequency 
^ 20 — 150 MHz ultrasound used for intravascular pho¬ 
toacoustic (IVPA) imaging Q . IVPA imaging is comple¬ 
mentary to traditional intravascular ultrasound (IVUS) 
imaging In IVUS imaging, ultrasound is used to 

interrogate the body, and an image is reconstructed from 
acoustic waves that are backscattered from blood vessel 
tissues. In IVPA imaging, a laser pulse is used to inter¬ 
rogate the body instead of ultrasound, and an image is 
reconstructed from broadband acoustic waves generated 
W light absorption events in various tissue components 
3 . Consequently, IVPA imaging is less affected by the 
light scattering in tissues, and therefore, can penetrate 
deeper as compared with purely optical imaging meth¬ 
ods [3-@. Moreover, the contrast of IVUS and IVPA 
imaging can be increased by exploiting harmonics of ul¬ 
trasound generated due to backscattering from nonlinear 
contrast agents Q. This fact additionally motivates our 
interest in the nonlinear acoustic interaction and optical 
detection of harmonics of ultrasound. 

Contrast in PA imaging arises from the natural varia¬ 
tion in the optical absorption of tissue components. The 
absorption cross-section of an optical antenna is many 
orders of magnitude higher than that of tissues. Conse¬ 
quently, recent attention has turned to the use of opti¬ 
cal antennae (plasmonic nanoparticles) as contrast agents 
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for IVPA imaging H^and all-optical sources of high- 
frequency ultrasound Q. 

Hereafter we theoretically demonstrate that in addi¬ 
tion to their ability to generate high-frequency ultra¬ 
sound, single optical antennae immersed into a bulk liq¬ 
uid can also be used to detect ultrasound. Ultrasonic 
transduction is via the detection of weak optical signals 
arising due to the Brillouin light scattering (BLS) from 
acoustic waves propagating in the liquid. This transduc¬ 
tion mechanism is different from that proposed in Ref. |^, 
where a large array of plasmonic metamaterials was used 
for non-resonant optical detection of ultrasound. 

Our idea is schematically illustrated in Fig. [T](a). The 
spectrum of light scattered from acoustic waves has a 
form of a triplet, consisting of the central Rayleigh peak 
and two Brillouin peaks shifted from the frequency of 
the incident light by the frequency of acoustic wave 
[Fig. mb)] m- In conventional approaches to detecting 
BLS, the Brillouin shift and the intensity of the Brillouin 
peaks are much smaller than the frequency and intensity 
of the incident light, respectively, and sophisticated in¬ 
strumentation |lll-ll^ has to be used to determine them 
experimentally. To address this problem, we numerically 
simulate the BLS effect from acoustic waves and show 
that the intensity of the BLS signal can be increased by 
exploiting higher-order plasmon modes supported by the 
dipole optical antenna. 

While there are many different optical antenna designs 
Qi , we focus our attention on dipole antennae that con¬ 
sist of a single nanorod made of silver. We are interested 
not only in the response of the dipole antenna to the 
fundamental resonance, but also to higher-order plasmon 
modes. Our focus on higher-order modes is to ultimately 
monitor backscattered acoustic radiation from biologi¬ 
cal tissue and bodily fluids, but also to detect nonlinear 
acoustic responses from such media [l^, which affords 
new sensing opportunities By using an ultrasmall 
optical antenna, rather than traditional mm-size, ultra¬ 
sound transducers, we see new opportunities for IVPA 
imaging Q. 
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FIG. 1. (a) Schematic of the detection of acoustic waves using 
an optical dipole antenna immersed into a liquid. An acoustic 
wave propagating in the liquid (blue background) is probed by 
an optical signal localised in the near-field-zone of the antenna 
(yellow rectangle). This signal originates from the nonlinear 
generation of the third harmonic of the frequency of the freely 
propagating light focused on the dipole optical antenna (de¬ 
picted by red arcs), (b) Typical BLS spectrum schematically 
plotted in a semi-log j/-axis. 


The existence of higher-order plasmon modes of dipole 
optical antennae have been theoretically predicted and 
experimentally observed by several groups (see, e.g., 
Refs. [iB-IH). These modes have been shown to be anal¬ 
ogous to the modes of a Fabry-Perot resonator. This 
analogy comes from the fact that these modes arise from 
multiple reflections between the two ends of the antenna 
[iM . Most significantly, lower radiative losses as¬ 

sociated with these modes lead to a higher local field 
confinement as compared with the fundamental mode. 
In particular, this property is attractive for the detection 
of weak dielectric permittivity fluctuations of the liquid 
around the antenna caused by ultrasound. We will ex¬ 
ploit this property in the remainder of the paper. 

Finally, although we are interested in linear and non¬ 
linear acoustic transduction, it is also important to stress 
that our antennae can demonstrate optical nonlinearities 
[3, [13 . By numerical simulations we demonstrate the 
possibility to optically excite the antenna by the near- 
infrared radiation but probe the acoustic waves with light 
at the frequency of the third harmonic (TH) of the op¬ 
tical excitation signal [Fig. HJa)]. This is achievable by 
designing the antenna such that the frequencies of one 
of its higher-order plasmon modes coincides with the fre¬ 
quency of the TH signal. Most significantly, a strong 
local field associated with the higher-order modes leads 
to a considerable enhancement of the TH signal in close 
vicinity to the antenna surface. This makes the antenna 
more sensitive to the fluctuations of the dielectric per¬ 
mittivity of the liquid that surrounds the antenna and 
transmits acoustic waves. Such nonlinear optical signal 
conversion can be advantageous in practical situations 


when the operation at shorter wavelengths is preferable 
due to specific absorption properties of biological tissues 
or liquids that are being insonified by ultrasound. More¬ 
over, the third harmonic generation (THG) itself can be 
used as an independent bioimaging modality [13 that 
can complement the IVPA imaging. 

The remainder of the paper is organised as follows. In 
subsection III Al we present a general theory of nonlinear 
acoustic wave propagation in bulk water, which is cho¬ 
sen as the model liquid surrounding the optical antenna. 
Based on this theory, we develop an acoustic finite- 
difference time-domain (a-FDTD) method, which allows 
calculating the spatial profiles of the scalar pressure field 
p{r,t) and the vector velocity field v{r,t) around the 
antenna. These profiles allow us to understand the prop¬ 
erties of the optical antenna to scatter acoustic waves. In 
subsection III Bl from the simulated profiles of the acous¬ 
tic fields we extract the acoustically-driven dynamics of 
the fluctuations of the dielectric permittivity of water, 
which is then used in our in-house optical finite-difference 
time-domain (o-FDTD) method modified to simulate the 
BLS spectra by taking into account the extremely slow 
dynamics of acoustics waves as compared with light. In 
subsection III Cl we demonstrate the enhancement of the 
BLS signal due to the higher-order plasmon modes of 
the antenna. Finally, in Section imi we discuss a poten¬ 
tial practical scheme of the excitation of the higher-order 
modes of the antenna by means of nonlinear optical ef¬ 
fects. We also present numerical details of the a-FDTD 
and o-FDTD methods in Appendix A and Appendix B, 
respectively. 


II. RESULTS AND DISCUSSION 

We explore the plasmon-assisted enhancement of the 
intensity of peaks arising in the spectral response of the 
optical antenna due to the BLS from acoustic waves, 
which propagate in bulk liquid that surrounds the an¬ 
tenna. To model the optical response of the antenna 
we employ the finite-difference time-domain (FDTD) 
method that is one of the standard numerical approaches 
employed in electrodynamics and photonics [^ . 

Finite-difference algorithms are also often used to solve 
acoustic wave propagation and scattering [ 23 . Concep¬ 
tually, the acoustic FDTD (a-FDTD) is similar to opti¬ 
cal FDTD (o-FDTD). This is because both approaches 
are based on first-order differential equations where the 
temporal derivative of one field is related to the spa¬ 
tial derivative of another field. Thus, it is theoreti¬ 
cally possible to merge the two FDTD algorithms to 
self-consistently solve the problem of BLS from acoustic 
waves in the presence of the optical antenna. 

In practice a rigorous self-consistent numerical solution 
that combines o-FDTD and a-FDTD is very challenging 
because the acoustic frequency (e.g., 10 MHz) and the 
speed of sound (e.g, 1500 m/s) are widely disparate from 
those of light (150—600 THz and 3x 10® m/s). Hence, the 
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numerical FDTD model effectively has two very different 
time scales. The first time scale is at the rapidly varying 
electromagnetic (light) frequency. The second time scale 
is the slowly varying acoustic frequency. As a result, the 
electromagnetic solution reaches a steady state solution 
before the liquid state has changed substantially due to 
the propagation of the acoustic wave. 

Inspired by the theoretical approach presented in 
Ref. I 27 I we circumvent the problem of different time 
scales by consideriM hypersonic acoustic waves that have 
GHz frequencies [l0|, instead of ultrasonic wave having 
several orders of magnitude lower frequencies. In many 
liquids, such waves can be excited by a laser pulse [281 
or piezoelectric film transducers . Importantly, in liq¬ 
uids the general properties of hypersonic waves are simi¬ 
lar to those of ultrasonic waves. Hence, the propagation 
of hypersonic waves can be studied by using conventional 
acoustic approaches 0- In the framework of our numer¬ 
ical model, the consideration of hypersonic waves instead 
of ultrasonic ones allows us to reduce, to some extent, the 
mismatch between the time scales of the o-FDTD and a- 
FDTD models. This eventually makes a full-vectorial 
three-dimensional 3D numerical model solvable when us¬ 
ing high-performance computers. However, a more sys¬ 
tematic analysis carried out in this paper has required us 
to additionally reduce the dimensionality of the acoustic 
model to two-dimensions (2D). 

For our purposes, we assume that our antenna is made 
of silver and it is surrounded by bulk water. This as¬ 
sumption is done without loss of generality because sil¬ 
ver is often used as the model material of optical anten¬ 
nae, and water is present in many practical situations 
including the propagation of ultrasound in biological ob¬ 
jects. Moreover, silver is used as the constituent material 
of nanoparticles employed as contrast agents for IVPA 
imaging Q- 

The propagation of acoustic waves in water (as well as 
in blood, fat, and other biological tissues and liquids) is 
governed by the equations of nonlinear acoustics derived 
from the full system of governing equations of fluid dy¬ 
namics [l^ ■ Indeed, differences between water and other 
fluids and materials is likely to be the ultimate sensing 
target of a final device and could, in principle, be treated 
with the approaches we outline below. 

The speed of the longitudinal acoustic wave cq in wa¬ 
ter (silver) is 1500 m/s (3650 m/s) and the density of 
water (silver) is 1000 kg/m^ (10500 kg/m^). These pa¬ 
rameters are taken into account by the model equations 
of nonlinear acoustics in liquids. 

However, silver also support transverse acoustic waves, 
the speed of which is smaller than that of longitudinal 
waves. The speed of transverse waves is a material pa¬ 
rameter taken into account in the equations of elastody- 
namics, which have to be solved if acoustic vibrations 
of the optical antenna are of interest. It is noteworthy 
that in optical antennae the frequency of these vibra¬ 
tions is ^ 10 — 20 GHz or even higher, which depends 
on the antenna dimensions [30l| . Hereafter, we neglect 


such vibrations, and thus we do not consider the elas¬ 
ticity equations in our model, which assumes that the 
water acoustic waves interact with the silver antenna as 
with a rigid body. This approximation is warranted be¬ 
cause the acoustic impedance of silver is much higher 
than that of water [^. Another reason to neglect the 
acoustic vibrations of the antenna is the fact that we con¬ 
sider GHz-range hypersonic waves instead of ultrasonic 
waves (see above). We project the results obtained for 
the hypersonic waves to the case of acoustic waves having 
lower frequencies, at which the acoustic vibrations of the 
ultrasmall antenna are negligible. This latter approxi¬ 
mation is also warranted because the dimensions of the 
antenna are small as compared with the wavelength of 
acoustic vibrations in the kHz and MHz spectral ranges. 
Indeed, acoustic vibrations at these frequencies are possi¬ 
ble in optomechanical micro-cavities whose characteristic 
dimensions reach 10 — 50 /rm [s^, i.e. ^ 30 — 150 times 
larger than an optical antenna. 


A. Model equations of nonlinear acoustics in 
liquids 

We begin with the continuity equation, which links the 
vector velocity vector v to the instantaneous density of 
liquid p as 


+ {v ■V)p + pV -v = 0, (1) 

where t is the time. 

The second equation is the momentum equation: 

p(^ + {v- V)») = 

-VP + rjV^v + (C + ^)V(V • v), 

where P is the thermodynamic pressure, 77 the shear vis¬ 
cosity, and C the bulk viscosity. 

The third, the thermodynamic equation of state, re¬ 
lates three physical quantities describing the thermody¬ 
namic behaviour of the liquid. These quantities can be 
P, p, and the temperature T. Alternatively, they can be 
P, p, and the specific entropy s. Normally, one applies a 
Taylor expansion about the equilibrium state, which to 
second order is 





( 3 ) 


where p, p', and s' are the dynamic pressure, density, and 
entropy, respectively. These quantities describe small dis¬ 
turbances relative to the uniform state of rest. The pa¬ 
rameters B jA and Cg are the nonlinear parameter of the 
medium and the ‘small signal’ speed of sound, respec¬ 
tively. 
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The fourth equation is the entropy equation 


pT (1^ + (u • V)s^ = + C(V • vf + 

^ 2 (4) 

p ^ _ 2 

2 v9a;j 5x1 3 '’■’5xk/ ’ 

where k is the thermal conductivity and is the delta 
function. The last term of the right hand side of Eq. |4]is 
written in Cartesian tensor notation. 

By keeping the first and second order terms in Eqs.[T]l4l 
we obtain the following nonlinear equations. The conti¬ 
nuity equation is 


dp' 1 dp^ 1 dC 

dt poCq dt Cq dt 

and the momentum equation is reexpressed as 


(5) 


dv „ 1 _ ‘iriy^dp „ „ 

Po— + Vp = 2 (C + 


dt 


PoCq 


dt 


( 6 ) 


where C = ^pou^ — is the second-order Lagrangian 
density. The entropy equation can be simplified to 


PoTo% = (7) 

dt 

where Tq the temperature of the medium at rest. The 
substitution of Eq. [3 to Eq. |3] and the use of thermody¬ 
namic relationships, one can write the equation of state 
as 


/ ^ _P_2_ (R\ 

^ cl pool 2pool Vcv Cp) dt' 

where Cv and Cp are the specific heat capacities at con¬ 
stant volume and pressure, respectively. We note that for 
a variety of liquids the constants used in Eqs. [T]l8] have 
been tabulated in the literature [s^- 

It is also important to derive the linearised acoustic 
equations. These equations serve as the departure point 
for the derivation of the discrete finite-difference ana¬ 
logues of nonlinear equations Eqs. [T]|8l In many cases 
it is instructive to start the analysis with the solution of 
the linear equations, and then to proceed with the solu¬ 
tion of the nonlinear equations, as the linearised equa¬ 
tions are easier to discretise and solve using the basic 
finite-difference formulations than their nonlinear coun¬ 
terparts [1^ [S^ . Additionally, the linear equations are 
also important because at the boundaries of the com¬ 
putational domain it is technically impossible to impose 
absorbing boundary conditions when using the nonlin¬ 
ear equations [2^. The absorbing boundary conditions 
is a standard approach that is used in the framework of 
the EDTD method to simulate the propagation of waves 
either in free space or in a medium in which the waves 
are fully attenuated before they can experience reflection 
from boundaries. 


Keeping only the term of the first order, the linearised 
continuity equation Eq. [T] becomes 

^ + V • (po^) = 0, (9) 

dt 

where po is the density of the undisturbed medium. Sim¬ 
ilarly, the momentum equation Eq. [5] becomes 

Po^ = -Vp- (10) 

The linearised equation of state Eq. |3] becomes 

P = cIp'. (II) 

Our discrete numerical model is based on the standard 
finite-difference stencils employed in the EDTD method 
[ 2 ^. The finite-difference analogues of both linear and 
nonlinear model equations are given in Appendix A. 


B. Fluctuations of the dielectric permittivity of 
liquid caused by acoustic waves 

The simulation of BLS from an acoustic wave requires 
one to define the time dependences of the fluctuations 
of the dielectric permittivity cased by the acoustic wave 
propagating in the liquid. This can be done by calcu¬ 
lating the profiles of the acoustic scalar pressure field p 
and the vector velocity field v. To start with, we cal¬ 
culate these profile for the linear case, i.e. we neglect 
the nonlinear acoustic effects. Once the behaviour of the 
dielectric permittivity in the linear case has been under¬ 
stood, we will take into account the nonlinear effects. 

So far, we have not discussed the particular design of 
the optical antenna. We assume that the antenna con¬ 
sists of a silver nanorod, which has a square-cross section 
w = 30 nm and length 340 nm. This choice of param¬ 
eters is due to optimisation using o-FDTD simulations, 
and is discussed in more detail in Section Ell 

In the a-FDTD model, the antenna is insonified by a 
plain quasi-monochromatic acoustic pressure wave prop¬ 
agating in bulk water. The peak amplitude of this wave is 
5 MPa, which is consistent with the pressure level used 
in nonlinear intravascular contrast imaging Q. In the 
time-domain, the profile of this wave is given by a very 
long Gaussian pulse modulated by a sinusoidal wave. The 
pulse length was 5000At where At is the time step in the 
a-FDTD algorithm. 

The left column of Fig. [3 shows the simulated profiles 
of |t;| in close proximity to the antenna for the acous¬ 
tic frequencies /a = I, 5, and 10 GHz. These profiles 
were obtained by Fourier-transforming the dependencies 
v{r,t) produced by the a-FDTD method. The contours 
of the antenna are outlined by the white 340 nm x 30 
nm rectangle, which also serves as the scalebar. We also 
calculate the displacement field which is the displace¬ 
ment of a particle from its equilibrium position due to 
the acoustic wave, at the same frequencies as above by 
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^ = f^v(r,t)dt (sTj. The calculation of ^ by means of 
the a-FDTD method is relatively simple when the inci¬ 
dent signal is quasi-monochromatic. 

As a next step, we calculate the divergence of the dis¬ 
placement field V • as shown in the right column of 
Fig. [21 The amplitude of dynamic fluctuations of the 
dielectric permittivity of water Se caused by the prop¬ 
agating acoustic wave is directly proportional to V • 
as shown in Ref. IH. This is because in the case of the 
plain longitudinal acoustic waves, the liquid is alternately 
compressed and stretched, which gives rise to alternating 
density changes. The knowledge of the changes in the di¬ 
electric permittivity is required to simulate the BLS from 
acoustic waves by means of the o-FDTD method. We es¬ 
timate that in our model the amplitude of Se is ^ 10“^. 

One can see that the spatial profile of V • ^ is quasi¬ 
uniform at the acoustic frequency /a = 1 GHz. This is 
because the wavelength of the acoustic wave in water is 
^ 5 times larger than the dimensions of the antenna. As 
a result, the acoustic wave is weakly scattered by the 
antenna. Consequently, the profile of the velocity v is 
quasi-uniform, which implies that in close proximity to 
the antenna V • ^ « 0. Hence, the profile of Se is quasi¬ 
uniform, too. The same result holds for lower frequencies 
(not shown). 

However, at higher frequencies the profile of V • ^ be¬ 
comes less uniform. This can be seen at ^ 5 GHz, i.e. 
when the acoustic wavelength 300 nm is close to the 
length of antenna. Moreover, at frequencies > 5 GHz 
the profile is no longer uniform (Fig. |21 right column). 
Consequently, if the frequency is < 5 GHz, which is the 
case of our study, we will abstract from the fact that we 
simulate the propagation of the GHz-range hypersonic 
waves. Instead, we will consider the normalised acoustic 
frequency /a//a.O) where /a,o is the frequency of the fun¬ 
damental harmonic of the incident wave. This approach 
allows us to project our results to the case of ultrasonic 
waves. 

We now consider nonlinear acoustic effects. To do so, 
in our a-FDTD software we turn on the acoustic nonlin¬ 
earity by setting /3 = 1 + ^ = 3.5, which corresponds 
to the nonlinear acoustic parameter of water [^. As 
shown in Appendix A, the nonlinearity is present in the 
ID model only. The input quasi-monochomatic pressure 
wave has the frequency /a,o and amplitude 5 MPa. The 
distortion due to the nonlinear propagation signal is de¬ 
tected at the output of the ID model and used as the 
excitation source in the same 2D model employed to pro¬ 
duce the results in Fig.|21 In general, the distorted signal 
carries energy at the fundamental frequency /a,o as well 
as at all harmonics generated due to the nonlinear in¬ 
teraction. It is noteworthy that the amplitude of the 
generated harmonics depends on the wave propagation 
distance in the liquid. In our model, this distance equals 
approximately one quarter of Xp = , where Aa,o 

is the wavelength of the fundamental wave and M the 
acoustic Mach number. The value of Xp corresponds to 
the distance at which the amplitude of the second non- 


pressure wave 



FIG. 2. Simulated spatial profiles of the velocity |t)| (left col¬ 
umn) and the divergence of the displacement field |^| for the 
acoustic frequency 1 GHz, 5 GHz, and 10 GHz. The white 
rectangle denotes the contours of the antenna and also serves 
as the scalebar with dimensions 340 x 30 nm. The black 
(yellow) colour denotes the zero (maximum) of the profile in¬ 
tensity. The direction of the incident monochromatic pressure 
wave is indicated by the arrow. Note that the velocity |t)| in 
front of (behind) the antenna increases (decreases), which also 
changes the profile of V ■ G This is because of the pressure 
doubling effect [3ll| . 


linear generated harmonics formally reaches 1/2 of the 
amplitude of the fundamental wave [3^ . 

The simulated spectrum of the linear acoustic wave 
is shown in Fig. |3Ka). It was obtained by Fourier- 
transforming the calculated time-domain dependencies of 
the pressure p{t) detected in the input port of the ID a- 
FDTD model. From Fig. EKa) and its inset, which shows 
the same spectrum plotted in the decibel scale, one can 
see only one peak at the frequency /a,o confirming the 
monochromatic nature of this wave. 

However, when the acoustic nonlinearity is taken into 
account, in Fig.jSKb) one can see new nonlinear generated 
waves appear in the spectrum at the frequencies 2/a,o 
(quadratic nonlinear effect), 3/a,o (cubic effect), etc. This 
spectrum was obtained by Fourier-transforming the val¬ 
ues of p{t) detected in the output port of the ID model. 
In addition to the harmonic waves, the nonlinear interac¬ 
tion leads to the creation of a dc component (i.e. zero fre¬ 
quency), which can be seen upon a close inspection of the 
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FIG. 3. Normalised acoustic energy spectrum of the incident 
quasi-monochormatic pressure wave detected (a) in the input 
port of the ID model and (b) in output port after the prop¬ 
agation in water. The frequency of the incident wave is /a,o 
and its amplitude is 5 MPa. The insets show the same spectra 
plotted in the decibel scale. Note that in the nonlinear spec¬ 
trum in total five new frequency harmonics can be seen above 
the —40 dB level. Furthermore, one can distinguish the dc 
component (zero frequency) generated due to the nonlinear 
acoustic effect. 

inset in Fig. [3Kb). We note that the relative amplitudes 
of the nonlinear generated waves are in good agreement 
with the predictions of analytical formulae presented in 
Ref. m. 

In the nonlinear case, we also extract the time-domain 
dependencies of the dielectric permittivity fluctuations 
Se{t) (not shown). In contrast to the linear case, these 
dependencies contain the contributions from the density 
changes caused by both the fundamental wave and its 
harmonics. The contributions of the harmonic waves are 
small as compared with those of the fundamental wave 
because their relative magnitudes follow the intensities 
of the spectral components shown in Fig. |3Kb). We use 
6e{t) in the simulations of light scattering from nonlin¬ 
ear acoustic waves. The results of these simulations are 
presented in subsection III Cl 


C. Plasmon-enhanced Brillouin light scattering 

We demonstrate a plasmon-enhanced light-acoustic 
wave interaction in the presence of the antenna. The 
a-FDTD simulations conducted in the previous Section 
III B I have produced the time dependences of the dielectric 


permittivity fluctuations at each discrete point inside the 
computational domain. We correlate these points with 
the discrete points of the computational domain in the 
o-FDTD method, as explained in Appendix B. 

We use the theory from Refs. [l|, [H, and [13 to es¬ 
timate the frequency of the fundamental (dipole) plas- 
monic mode of the antenna: ^ 2100 nm. The theory also 
predicts that this antenna should support two higher- 
order plasmonic modes at ~ 700 and ^ 500 nm, respec¬ 
tively. 

Hereafter, we also assume that the acoustic waves in 
water are probed at the frequency of the first higher- 
order plasmonic mode at ~ 700 nm. The reason for this 
is twofold. As will be shown in Section imi the electric 
field profiles of the higher-order plasmonic modes exhibit 
Fabry-Perot-like behaviour and are tightly confined to 
the surface of the antenna. These modes bounce back 
and forth as they experience multiple reflections from the 
ends of the antenna. This behaviour is favourable for the 
detection of acoustic waves because the multiple passing 
of light along the antenna surface effectively increases 
the interaction time of light with the acoustic wave. Fur¬ 
thermore, the tight confinement of the electric field gives 
rise to an increased sensitivity of the plasmon resonance 
to any change in the dielectric permittivity occurring in 
close proximity to the surface of the antenna [3^ . 

It is important to mention that for numerical tractabil- 
ity we have to make several simplifications, which do not 
interfere with our analysis but drastically decrease sim¬ 
ulation time and simplify the analysis of the plasmon- 
enhanced BLS effect. 

Firstly, we limit ourselves to the 2D o-FDTD model. 
Apart from the fact that the a-FDTD model is also 2D, 
this simplification is warranted because 2D o-FDTD sim¬ 
ulations produce qualitatively the same linear spectra as 
those produced using a 3D model (see, e.g., Ref. 113). 

Secondly, we assume that the optical signal exciting the 
antenna does not generate additional acoustic waves in 
water. This is because an efficient generation of acoustic 
waves in water is possible only when one fulfils certain 
conditions 0 , 113 . However, those conditions are not met 
in our model because their fulfilment is out of the scope of 
this work. Furthermore, we also neglect the impact of the 
radiation pressure, electrostriction, thermal expansion, 
and other effects that potentially can give rise to acoustic 
vibrations of the antenna due to the interaction with the 
optical radiation. 

To start with, we simulate the effect of BLS from 
a quasi-monochromatic sinusoidal acoustic wave of fre¬ 
quency /a,o propagating in the linear medium. In the 
a-FDTD software, we turn off the nonlinearity by setting 
/3 = 0. We run the a-FDTD simulation in parallel with 
the o-FDTD simulation, where the latter receives time- 
dependent information about the dielectric permittivity 
fluctuations of water caused by the propagating acoustic 
wave. Two batches of simulations are conducted: one 
with the optical antenna in bulk water and another one 
in bulk water without the antenna. The results are pre- 
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sented in Fig.|4lja). In the case without the antenna (blue 
curve) one can see a typical BLS spectrum with the the 
central Rayleigh peak and two Brillouin peaks shifted 
from the frequency of the incident light by the frequency 
of acoustic wave /a,o- 

Importantly, in the presence of the antenna one can 
see that the intensity of Brillouin peaks is enhanced (red 
curve). This ^ 35-fold increase with respect to the BLS 
intensity without the antenna is attributed to the plas- 
mon resonance of the antenna. To verify this, we con¬ 
ducted an extra simulation (results not shown here) in 
which the frequency of the optical signal that excites the 
antenna was off-resonance. As expected, that simula¬ 
tion produced virtually the same result as without the 
antenna. 

We also note the existence of two additional Brillouin 
peaks shifted from the frequency of the incident light 
by 2/a,o [Fig. lUJa)]. As the acoustic nonlinearities were 
neglected in the simulation, the appearance of these very 
weak peaks is attributed to the cascaded scattering of the 
optical signals at the frequency ±/a,o from the acoustic 
wave. In this case, the optical waves at ±/a,o act as the 
incident light but the signals at ±2/a,o are the Brillouin- 
shifted waves, which are very weak and visible only in 
the presence of the antenna. 

Finally, we simulate the BLS from nonlinear acoustic 
waves in the presence of the antenna. We repeat the 
same self-consistent acoustic and optical simulations as 
above, but this time we take into account the nonlinear 
acoustic properties of water. The result of this simulation 
is presented in Fig. |ll[b). In contrast to the BLS effect 
from the linear wave [Fig. HJa), red curve], in Fig. SKb) 
one can see multiple Brillouin peaks corresponding to 
the scattering from nonlinear generated harmonics of the 
incident acoustic wave. Importantly, the enhancement 
of the BLS intensity by the antenna occurs only at the 
frequencies of the Brillouin peaks and it does not affect 
a broad Rayleigh-wing background. This can be seen by 
inspecting the difference between the red and blue curves 
in Fig.|4Kb) that correspond to the cases with and without 
the antenna, respectively. This selective enhancement 
will help to resolve the fine structure in BLS spectra. 

It is noteworthy that the optical spectrum in Fig. SKb) 
resembles a Brillouin frequency comb similar to that in 
Ref. 1^ However, the physical origin of this frequency 
comb is different. This comb is generated due to the 
BLS from the nonlinear acoustic waves as well as the 
subsequent plasmon-assisted enhancement of the BLS 
peaks. We verihed that the intensity of the individual 
optical spectral components in Fig. |4][b) corresponds to 
the intensity of the harmonics in the acoustic spectrum 
in Fig. [3Kb). We note that the spectral composition of 
the frequency comb in Fig. SKb) can be controlled by tai¬ 
loring the nonlinear acoustic interaction in water. This 
functionality can be exploited in recently proposed Bril¬ 
louin cavity opto-mechanic systems combined with mi¬ 
crofluidic devices (40l| . 



FIG. 4. Simulated plasmon-enhanced intensity of light scat¬ 
tered from (a) the monochromatic sinusoidal acoustic wave 
and (b) nonlinear generated acoustic waves. Note a 35-fold 
enhancement due to the presence of the nanorod in Panel (a). 


III. DESIGN OF THE OPTICAL ANTENNA 


We discuss the possible design of the optical antenna 
and one of the possible practical schemes for its optical 
excitation. In the discussion above, we have already men¬ 
tioned the existence of higher-order plasmon modes sup¬ 
ported by the antenna. We also assumed that in our sim¬ 
ulations of the plasmon-enhanced BLS effect, the acoustic 
waves were probed by light whose frequency corresponds 
to the resonance frequency of the excited higher-order 
mode. The frequency of this mode was estimated using 
simple effective wavelength theory. 

Recall that the higher-order modes of the antenna are 
analogous to Fabry-Perot modes because they arise from 
multiple reflections between the two ends of the silver 
nanorod. It is noteworthy that different Fabry-Perot 
structures and other types of optical resonators are of¬ 
ten used in many practical situations involving the de¬ 
tection of acoustic waves with light (see, e.g., Ref. |4l| 
and references therein). This is because the presence of 
the resonator effectively increases the interaction time of 
light with the acoustic waves. We also note that the reso¬ 
nant excitation of plasmons has been used to additionally 
improve the resolution and sensitivity of the BLS spec¬ 
troscopy and relevant techniques (see, e.g., Refs. 113 and 
113). Consequently, the higher-order modes effectively 
combine both the resonant excitation of plasmons with 
the Fabry-Perot-like behaviour. 

We also remind that the far-held optical radiation 
properties of the higher-order plasmonic modes of the 
antenna are poor as compared with those of the dipole 













modes. By virtue of the principle of reciprocity [l|, this 
implies that the efficiency of the excitation of the higher- 
order modes from the far-held region is also relatively 
low. Consequently, these modes typically have to be ex¬ 
cited by near-held optical radiation, for example, by us¬ 
ing aj)oint-like quantum nano-emitter such as a quantum 
dot [l|. 

However, there are alternatives to the near-held exci¬ 
tation scheme involving nano-emitters. We demonstrate 
that the higher-order modes can be excited with the far- 
held radiation via nonlinear third harmonic generation. 
Previously, Abb et. al. [1^ demonstrated a nonlinear 
control of the higher-order modes in asymmetric nanorod 
antenna with a gap. Here, we show that the light inten¬ 
sity enhancement in the near-held region associated with 
the excitation of these modes leads to a considerable in¬ 
crease of the THG effect as compared with the case with¬ 
out the antenna. This distinguishes our approach from 
that used in many previous works that exploit dipole 
plasmonic modes to enhance nonlinear optical effects at 
the nanoscale (see, e.g., Refs. Idsl - ldsll . 

We continue considering the silver antenna with a 
square cross-section re = 30 nm and length L = 340 
nm. However, since optical simulations in the absence of 
acoustic waves are signihcantly faster, we conduct rigor¬ 
ous 3D numerical simulations of the linear and nonlinear 
response of the antenna. The conhguration of our model 
is shown in Fig. [SJa). The antenna is embedded into a 
uniform loss-less dielectric medium with an effective di¬ 
electric constant Ceff = 1-75. Whereas this value is close 
to the optical dielectric permittivity of water, it also ac¬ 
counts for the presence of a dielectric substrate of the 
antenna or an optical fibre use to guide light towards the 
antenna. These components can be present in a realistic 
device. 

Figure [SKb) (red curve) shows the simulated linear 
spectrum of the antenna. One can see a broad dipole 
mode peak at ~ 2000 nm and two narrow higher-order 
mode peaks at 750 nm and 550 nm. Remarkably, the 
spectral position of the observed peaks is in good agree¬ 
ment with the approximate resonant wavelengths pro¬ 
duced by a simple theoretical model based on the effec¬ 
tive wavelength. 

In Fig. [S](c) we present the profiles of the normalised 
electric field intensity through the centre of the cross- 
section, i.e. along the j/-axis in Fig. [SKa). The red solid 
curve denotes the profile for the fundamental mode. The 
blue dashed (black dash-dotted) curve denotes the profile 
of the higher-order mode at 750 nm (550 nm). The rect¬ 
angle above Fig. [5Kc) schematically indicates the length 
of the antenna. One can see that the electric field profile 
of the higher-order plasmon modes resembles that for a 
Fabry-Perot resonator. 

Now we study the nonlinear optical THG effect arising 
due the optical antenna. We conduct simulations of the 
THG effect by using a 3D nonlinear FDTD method ; 
for more details see Appendix B. We neglect the nonlin¬ 
ear susceptibility of silver and for the surrounding ef- 




FIG. 5. (a) Illustration of the numerical simulation config¬ 
uration for the investigation of the THG in the antenna. A 
plane electromagnetic wave is normally incident from the left 
on the antenna. The electric field vector of the plane wave is 
parallel to the long axis of the antenna. The generated TH 
signal, as well as the signal and the fundamental frequency, 
are probed from the right in a point located 50 nm off the 
middle point of the nanorod. (b) Left j/-axis: Normalised to 
unity linear optical intensity spectrum of the antenna. Right 
j;-axis: THG enhancement spectrum, (c) Normalised simu¬ 
lated electric field intensity through the centre of the antenna 
cross-section. The rectangle above Panel (c) schematically 
indicates the length of the antenna. 

fective dielectric medium we assume that = 1x10“^^ 
m^/V^, which is close to the realistic value for water and 
other liquids [1^. Although the nonlinear susceptibility 
of silver can be higher than that of water, its neglect is 
justified by the fact that silver may be oxidised in the 
presence of water, which is expected to significantly re¬ 
duce its contribution to the nonlinear response. 

The antenna is excited in the 1200 — 2500 nm spectral 
range with the amplitude of the electric field of Ey = 
1.0 X 10® V/m corresponding to the optical intensity of 
^ 2.3 GW/cm^. In our simulations, the intensity of the 
incident light is constant at all frequencies. This is a 
convenient theoretical approach that allows simplifying 
the analysis of results as compared with the case of the 
more realistic frequency-dependent intensity. 
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We study the THG enhancement defined as a ratio of 
the intensity of the generated TH signal in the presence 
of the antenna to that without it. We expect to observe 
a very large enhancement of several orders of magnitude 
because the intensity of the TH signal simulated in the 
absence of the antenna should be negligibly small due 
to a very short nonlinear optical interaction length. In 
our model this length equals 500 nm. Indeed, as shown 
by open blue circles in Fig. [5jb), the THG enhancement 
reaches a remarkably large value ^ 10 ® at the frequency 
of the first higher-order mode located at 750 nm. This 
also corresponds to the conversion efficiency 77 = 0.033%, 
which was calculated as the I^ui/Iui being the inten¬ 
sity of the TH signal and lu) the intensity of the dipole 
mode of the optical antenna. We also note that even 
though the wavelength of the second higher-order peak 
(550 nm) is detuned from the wavelength of the TH sig¬ 
nal, the THG enhancement spectrum at 550 nm remains 
noticeable. 


IV. CONCLUSIONS 


that many high energy sonochemical reactions achievable 
with ultrasound [s^ exploit the effect of cavitation, i.e. 
the formation, growth and collapse of bubbles in a liquid 
Remarkably, the nonlinear acoustic coefficient of a 
bubbly liquid can be many orders of magnitude higher 
than that of the same liquid but without bubbles. Al¬ 
though the consideration of such strong acoustic nonlin¬ 
earities goes far beyond the scope of this paper, we en¬ 
vision that the optical antenna can be used to facilitate 
the detection of highly nonlinear acoustic waves. 
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We have theoretically demonstrated the possibility to 
employ dipole optical antennae to enhance the scatter¬ 
ing of light with nonlinear acoustic waves propagating in 
liquids. Although one can use the same optical antennae 
as those already employed in plasmon-enhanced photoa¬ 
coustic imaging of biological tissues, we have showcased 
a different antenna design - single 340 nm-long silver 
nanorods that are long enough to support both funda¬ 
mental and higher-order plasmon modes. We have shown 
that the higher-order plasmon modes of these nanorods 
lead to a strong interaction of light with acoustic waves. 
This effect originates from both the near-field enhance¬ 
ment due to the plasmon resonance and the Fabry-Perot- 
like behaviour of the antenna. Furthermore, we have 
shown the possibility to excite the optical antenna at one 
optical frequency but detect the acoustic waves at an¬ 
other optical frequency. This functionality is achievable 
by means of the nonlinear optical Kerr effect exhibited by 
water, which was considered as an ideal liquid that sur¬ 
rounds the antenna and also transmits the acoustic wave. 
We envision that this optical excitation scheme can find 
practical applications in biomedical imaging where dif¬ 
ferent tissues and liquids have widely disparate optical 
absorption spectra. We also presented a hybrid acousto- 
optical finite-difference time-domain algorithm that al¬ 
lows one to investigate the interaction of light with non¬ 
linear acoustic waves at the nanoscale. 

Furthermore, the optical antennae can be used to addi¬ 
tionally increase the previously demonstrated benefits of 
the BLS spectroscopy technique from plasmon resonances 
in metallic nanostructures [^, |^, . Ultrasound is also 

an important tool in food technology (see, e.g., Ref. H^, 
and optical antennae may be useful as a component of an 
all-optical system for noninvasive quality control of dif¬ 
ferent processes involving ultrasonic waves. We also note 


Appendix A: The acoustic FDTD model 

The derivation of the discrete nonlinear acoustic equa¬ 
tions is difficult as compared with the case of the lin¬ 
earised equations. Consequently, it is instructive to de¬ 
rive the discrete analogues of the linearised equations 
first, and then to use the resulting equations as the de¬ 
parture point for the derivation of the nonlinear equa¬ 
tion. Using the FDTD notation proposed by Yee [ 2 ^ , in 
the two-dimensional (2D) x — y space our linear FDTD 
equations for the vector velocity v are written as 

^ b”(bj) -p”(bj -1)] 

^ K(bj) -p”(*- 1,7)] 

(Al) 

where the index n denotes the current iteration number. 
For the scalar pressure field p we obtain 




-f 


Ay 

Kp 


+ 1 ) - 


(A2) 


Eqs. IA1IIA2I include the co-ordinate dependent coef¬ 
ficients Kv = At/po and Kp = AtpoCg, where At = 
Ax/(-\/2co,max) Is the time step calculated assuming that 
Ax = Ay and that co,max is the maximum speed of sound 
propagation in the liquid. Eqs. IA1IIA2I are solved it¬ 
eratively in a fashion similar to electromagnetic FDTD 
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solvers [1^ . The information about how to model acous¬ 
tically rigid bodies can be found in Ref. [H. 

It is noteworthy that in our simulations the nonlin¬ 
ear acoustic effects are taken into account as a one¬ 
dimensional (ID) model. This is because in many practi¬ 
cal cases it suffices to solve a ID wave propagation prob¬ 
lem in order to explain experimental results |5l|. In addi¬ 
tion, we do not expect that the presence of an ultra-small 
antenna to significantly perturb the strength of the non¬ 
linear acoustic effects in water. Consequently, we use the 
ID nonlinear model to produce the frequencies and am¬ 
plitudes of the harmonics of the original acoustic wave. 
As a next step of our simulations, these data are used in 
the linear 2D model. 

Following Ref. by using the first-order differencing 
in time and space, we derive the following finite-difference 
equations for the nonlinear ID case 

(k) = v^-^{k) + 5b [p{kT -P"{k + 1)] 

+ 1) - - p^-\k + I) +p^-\k)\ 

(A3) 

where the constants cq, po, Ci p are co-ordinate de¬ 
pendent, 5b = AtAz/po and we also chose another co¬ 
ordinate index variable, k, and another step in space, 
Az, in order to avoid the possible confusion with their 
counterparts in the 2D model above. Importantly, as 
compared with the linear algorithm, the values of pres¬ 
sure p at the two time steps - and - must 

be additionally retained in the computer memory. These 
extra values are required to update Eq. IA3I as well as 
Eq. IA4I which reads 


p“+i(fc) = - - {[coAzp""^(A:) - 2coAzp"(fc)] poC 

-f c^A^p^ (fc - I) - (k) 

+ CpAt Azp(fc)po 

- ;9AtAzp"(fc)p''-i(fc) + /3AtAzp"(fc)2} 


(A4) 


where 6 = [cgA^poC - cgAtA^po - +/3AtA*p“(fe)], 

/3 = l+ SandC'=^('A_J_y 
2A POC* \^Cv Cp y 

The following material constants of water were used in 
our simulations: the shear viscosity p = 1.002 x 10“^ Pa s, 
bulk viscosity ( = 3.09 x 10“^ Pa s, specific heat capac¬ 
ity at constant pressure Cp = 4.182 x 10^ J/(kg K) and 
constant volume Cv = Cn/1.33, and thermal conductivity 
AC = 0.597 W/(m K) [13. 

Eqs. IA31IA4I are cumbersome but their solution by 
means of the standard FDTD procedure is straightfor¬ 
ward and it does not suffer from long-term instabilities 
possible in some nonlinear FDTD schemes [13 . However, 
to achieve a better stability one needs to use a more strin¬ 
gent time step, as directed by the Courant stability fac¬ 
tor At = 0.5Az/co,max- Naturally, the standard Courant 


factor (without the 0.5 pre-factor) works very well in the 
case of linearised equations, i.e. when /3 = 0. 

The selection of the right value of Az is essential 
for conducting a high-accuracy simulation. (We assume 
that in the 2D model the steps in space along the x— 
and y— directions are both equal to Az.) Normally, 
a high accuracy is achievable by choosing Az < A/10 
where A is the shortest acoustic wavelength present in the 
model. For instance, for the acoustic frequency 1 GHz 
and Az = 2 nm one obtains 750 spatial points per wave¬ 
length. Whereas a 50 times larger step in space would 
already suffice to ensure a high accuracy, Az = 2 nm is 
used in this paper because we also model the antenna 
whose cross-section equals 30 nm, which is equivalent to 
15 spatial points. (Simulation of the antenna would be 
inaccurate if the cross-section were ~ 10 spatial points). 

Furthermore, we note that in the acoustic simulations 
the Courant stability criterion produces At ~ 5 x 10“^^ 
s when Az = 2 nm. It is also well-known that at least 
3 — 4 oscillation periods of the sinusoidal acoustic wave 
have to be simulated to obtain the correct result. To 
meet this condition, for the acoustic frequency 1 GHz 
one needs to conduct ~ 10^ iterations, which is equiva¬ 
lent to ^ 15 minutes of CPU time when the simulation 
is run on a multi-CPU workstation. Now let us assume 
that the acoustic frequency equals 1 MHz. As we still use 
Az = 2 nm required to correctly model the antenna, we 
obtain 750000 spatial points per wavelength and we need 
to conduct ^ 10^ iterations. This requires ^ 2500 hours 
of CPU time. Such a long simulation time is prohibitive 
even though a significant acceleration might be achieved 
by running the simulation on a supercomputer. This 
discussion additionally justifies our approach to model 
GHz-range hypersonic waves and project the results to 
the case of ultrasonic waves having lower frequencies. 


Appendix B: Connecting the acoustic and the 
optical FDTD models 

In this section we show how to connect the a-FDTD 
model with the o-FDTD model. We also discuss the non¬ 
linear o-FDTD used to simulate the purely optical prop¬ 
erties of the antenna. We start with the basic equations 
of the o-FDTD model: 

_„2I^=VxE(r,i), (Bl) 

at 


and 


+ J{r,t) = V X H{r,t), (B2) 

at 

where E{r,t) is the electric field, H{r,t) is the magnetic 
field, and J{r,t) is the electric current density required 
to implement the optical Kerr nonlinearity |46l |. The di¬ 
electric permittivity term e also enters this equation and, 
as the optical properties of silver are modelled by using 
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a Drude model [1^, in the absence of acoustic waves it 
can be assumed to be independent of r. 

However, in the presence of acoustic waves in agree¬ 
ment with Ref. l35| the total permittivity etot(j’,t) = 
e -I- 5f.{r, t) has to be introduced. The fluctuations of the 
dielectric permittivity due to the propagating acoustic 
waves 6e{r,t) is calculated from the values of the vector 
velocity field v{r,t) produced by the a-FDTD method. 

Consequently, Eq. IB2I reads 

[e + Se{r,t)] = V x H{r,t)- (B3) 

We remind that in the simulations of the plasmon- 
enhanced BLS we assume that outside the antenna 5e 
is independent of the space coordinate r. However, in 
Eq. IB3I we keep this dependence for the sake of general¬ 
ity. 

The presence of 6e{r,t) in Eq. IB3I is the only princi¬ 
pal difference between the o-FDTD method used in this 
work to simulate the BLS effect and the standard electro¬ 
magnetic FDTD method [ 2 ^ , which also used in Section 
m to investigate the purely optical properties of the an¬ 
tenna. Consequently, it is straightforward to derive the 
discrete analogues of Eqs. IBlI and IB3I 

The iterative solution of the discrete equations is only 
possible by properly addressing the mismatch between 
the frequency of light and the frequency of the acoustic 
waves. As has been shown in the main text, the dielectric 
permittivity fluctuations are quasi-uniform at frequencies 
below or equal to 1 GHz. This observation drastically 
simplifies the stimulations. We additionally accelerate 
the simulations by taking advantage of a slow change of 
the dielectric permittivity as compared with the optical 
response of the antenna. We run the o-FDTD simulator 
at discrete values of time at which the dielectric per¬ 


mittivity change is large enough to modify the optical 
response. 

The consideration of the nonlinear optical Kerr effect 
is also not straightforward because nonlinear effect are 
not taken into account in the standard o-FDTD formu¬ 
lations [ 2 ^ . As shown in Ref. |4^ in the time domain the 
polarisation associated with the Kerr nonlinearity can be 
written as 

PKerrit) = eox‘'^^\Eit)\^E{t). (B4) 

We use one of the noniterative updating schemes pre¬ 
sented in Ref . 1^ . To fulfil the causality principle of non¬ 
linear optics , we discretise Eq. IB4I as 

PKerr(tn+l) = \E{t^)\^ E^U+i), (B5) 

where n stands for the time in step used in the o-FDTD 
model and is the third-order nonlinear susceptibil¬ 
ity of the medium. Then we substitute Eq. IB5I into the 
constitutive equation D{r,t) = €eoE{r,t) -|- PTie„{r,t), 
where we neglected the dependence of the dielectric per¬ 
mittivity e on the fluctuations caused by the acoustic 
waves because these fluctuations are irrelevant to the 
study of the nonlinear optical properties of the antenna. 
Finally, we obtain the following noniterative electric-field 
updating equation 

= D"+V (eeo + , (B6) 

where q = {x, y, z) denotes the components of the electric 
field vector E in the Cartesian coordinate system and the 
discrete time tn has been replaced by the index n. 
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